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predicts pinning of the contact line at rectangular ridges perpendicular to flow for contact an- 
gles e > 45°. While for 6 = 30°, 6= 40° (no flow) and d = 60° (flow) all methods are found 
to produce data consistent with the CF criterion, at = 50° the numerical experiments provide 
different results. Whilst pinning of the liquid front is observed both in the LB and CFD simu- 
lations, MD simulations show that molecular fluctuations allow front propagation even above 
the critical value predicted by the deterministic CF criterion, thereby introducing a sensitivity 
to the obstacle height. 

Introduction 

In the recent years, with the rapid technological progress in the production of micro- and nano- 
channels, the understanding of fluid flow on the nanoscale^i^ has become crucial for modern nan- 
otechnology (such as the "lab on a chip" and related microfluidic devices) as well as for various 
applications of porous materials, flow in biomembranes, etc. A key problem is the description 
of fluid flow in narrow channels with wettable walls. Such channels are ubiquitous in cells and 
living matter but have been also successfully produced from synthetic materials in recent years .- 
Thus, planar nanochannels, fabricated by silicon-based technology, can provide an attractive con- 
figuration for fundamental studies like filling kinetics,^ hydrodynamics in confinement, and for 
molecular separation processes in biology.- Indeed, besides practical applications, microfluidics 
also raises a challenge to basic research since the continuum description of fluid flow goes under 
question whenever discreteness of matter comes into play, that is, at length scales comparable to 
the molecular size. 

To gain insight into the transport mechanisms involved in fluid flows, many researchers have 
studied the problem using a variety of computer simulation methods, and most notably Compu- 
tational Fluid Dynamics (CFD), Lattice-Boltzmann Equations (LBE), and Molecular Dynamics 
(MD) methods. The classical continuum theory based on Navier-Stokes (NS) equations assumes 
that state variables do not vary appreciably on a length scale comparable to the molecular free 
path. This conjecture has been challenged by both experiment^ and the earliest MD computer 
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simulations- which indicate the existence of significant density fluctuations of the liquid normal 
to the solid wall. Against expectations, however, some MD studies^ on PoiseuUe flow demon- 
strated that the classical NS description may be used for modeling capillary flow in channels with 
diameter of several molecular sizes and greater, while this has been challenged by another study.— 
To these controversial results one should add data, produced by the LBE method— which has 
gained increased prominence in recent years due to its efficiency and proximity to the basic as- 
sumptions of the NS constitutive equations. Results from the LBE approach indicate that there 
are several microfluidic situations, in which molecular details, although non-negligible, can still 
be given a meso-scopic, rather than fully atomistic, representation, without affecting the basic 
physics. i2ii3,i4,i5 ti6J2Jl On general grounds, this can be expected to be the case whenever molecu- 
lar fluctuations do not play any major role. However, as the physics of micro/nanoflows progresses 
towards increasingly demanding standards, qualitative expectations need to be complemented and 
possibly tested against quantitative assessments. 

In view of the diversity of methods and the plethora of controversial results from the computa- 
tional modeling of flow in the sub-micrometer range, a test of the adequacy and reliability of the 
principal approaches is highly warranted. In the present investigation we perform such a test by 
comparing the results from three of the most broadly used methods - CFD, LBE, and MD - focus- 
ing on a generic case, the capillary filling of a wettable narrow channel by spontaneous imbibition 
of a simple fluid. 

The aim of this comparative study is as follows: (i) we test the reliability of the three simulation 
methods with respect to the capillary filling of a wettable nanochannel, by comparing the results 
against the Lucas - Washburn (LW) theoretical prediction; (ii) we analyse the effect of the presence 
of an isolated corrugation on the flow dynamics and test the numerical results against the theoretical 
Concus-Finn criterion; (iii) the direct comparison between fluctuating (MD) and non-fluctuating 
(CFD and LB) approaches, permits to highlight the role of thermal fluctuations on the universality 
of the CF criterion. 



3 



S. Chibbaro et al. Capillary filling with wall corrugations 

Theoretical background 

In a horizontal capillary under steady state conditions (i.e., in the absence of gravity and transient 
inertial effects), the capillary imbibition pressure is balanced by the viscous drag of the liquid. 
Simple analysis of this process leads to the Lucas - Washburn equation (LW)— which relates the 
distance of penetration z{t) of the fluid front at time t to the capillary size H, the viscosity rj and 
surface tension 7 of the liquid, and the static contact angle 6 between the liquid and the capillary 
wall. If slip velocity and deviations from Poiseille flow are neglected, in the late asymptotic regime 
the distance z{t) travelled by the moving interface in the channel (having zo as initial coordinate) 
is given by 

z{t)^-zl= ^ ^ ^ ' ct, (1) 

where H denotes the height of the channel and the constant c =1/3 for a slit-like capillary. Equa- 
tion[I]can be recast in dimensionless form as z = z/H and i = t/ tcap, where tcap = Hrj / 7, becoming: 

z{i)^-zo = cos{e)ci. (2) 



It is apparent that in these dimensionless units, time-penetration depends only on the value of the 
contact angle 6, regardless of other fluid parameters and geometrical details. This facilitates the 
comparison between different simulation methods. 




H 



Figure 1 : Scheme of the microchannel geometry 



In this comparative study we explore capillary filling in the presence of topological obstacles 



(rectangular ridges) on the channel wall (Figure 1 ). Since capillary filling is mainly determined by 
the three-phase boundary ("contact line") between liquid, wall and vapor, the motion of the contact 
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line offers a much more stringent test of the various simulation methods and offers a possibility to 
assess their shortcomings and advantages. A key problem in this regard is the pinning of the contact 
line, due to a geometric singularity in the meniscus stability, like the presence of obstacles whose 
sides make an angle 2 a with respect to the wall, broadly known as Concus-Finn condition.-^i*2^i^ 
According to the Concus-Finn criterion, there is no filling due to meniscus instability, provided the 
contact angle fulfills the following condition: 



For a rectangular obstacle, this means that a contact line will be pinned at its trailing edge, once > 
45^, regardless of the obstacle height. While this condition follows from the detailed mathematical 
analysis— of the liquid surface stability, in fact it goes back to the thermodynamic foundations of 
wetting, as established in the early works of Gibbs.^^ The possibility for capillary driven flow is of 
major importance for numerous applications, e.g., modem fuel cells,— therefore it is not surprising 
that this criterion has been tested even in space experiments on the Shuttle-Mir complex. 

Modelling microfluids: from continuum mechanics to molecular 
dynamics 

Currently, intensive efforts to get deeper insight and understanding of flow in microchannels are 
carried out by researchers using a variety of computer modeling approaches. In principle, con- 
tinuum fluid mechanics provides the most economical description of micro flows. However, this 
approach fails to describe a series of phenomena, such as near-wall slip-flow, in which the contin- 
uum assumption goes under question due to the onset of molecular effects, especially near solid 
walls. Molecular Dynamics provides a high degree of physical fidelity, at the price, however, of 
a very susbtantial computational burden. The lattice kinetic approach is well positioned to offer a 
good compromise between the physical realism of molecular dynamics and the computational effi- 




- a. 



(3) 
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ciency of continuum mechanics. -^^^ The lattice kinetic approach is finding increasing applications 
to microfiuidic problems, as it permits to handle fluid- wall interactions at a more microscopic level 
than the Navier-Stokes equations, while simultaneously reaching up to much larger scales than 
Molecular-Dynamics.— ' ^'^ ' ^ However, the LB approach is subject to a number of limitations, 
such as the existence of spurious currents near curved interfaces, as well as enhanced evapora- 
tion/condensation effects due to the finite-width of the interface.-^ Ideally, one would combine the 
three methods within a synergistic multiscale procedure, using each of them wherever/whenever 
appropriate, in different regions of the microflow. Various two-level (continuum-MD, LB-MD and 
continuum-LB), are already available in the literature, while, to the best of our knowledge, three- 
level coupling are just beginning to appear.— A detailed and comparative assessment of merits and 
downsides of the three levels of description, is therefore of great importance to the ultimate goal 
of developing efficient and robust multiscale coupling methods for complex microfiuidic fiows. In 
this work, we shall present a comparative study of capillary filling in the presence of wall corru- 
gations, using all of the three methods mentioned above, namely a finite-volume CFD software 
package, a LB model with non-ideal fluid interactions, and a MD computer program. Since all of 
these methods are by now standard, in the following we shall provide just a cursory description, 
leaving the details to the vaste literature on the subject. 

Computational Fluid Dynamics 

Our CFD approach is based on the CFD-ACE-i- software package (release 2008), as a multiphysics 
commercial tool including geometry modeling, grid generation and results visualization.-^^ The 
spontaneous capillary filling in micro/nano channels is reproduced by means of the VOF (Volume 
of Fluid) scheme, based on the hypotheses of incompressibility, no-interpenetration and negligible 
localized relative slip of the two fluids in contact. In order to describe the transport of the volume 
fraction ^ of one of the two fluids in each cell (0 thus ranging from to 1), the Navier-Stokes 
equations for the fluid velocity are augmented with a scalar transport equation for the volume 
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fraction: 



50 

'di 



+ V-(m0) =0 



(4) 



where t is time, V is the standard spatial gradient operator, and u is the velocity vector of the 
fluid. The composition of the two fluids in the mixture determines for each computational cell the 
averaged value of physical properties such as density and viscosity. Any volume- specific quantity 
is evaluated in accordance with the following expression: 



where Cb is the volume-averaged quantity, (D/ (resp. cOg) is the value of the property for one liquid 
(resp. gas). However, being the averaged volume fraction of fluid in each cell, the definition of 
the interface between the two fluids in that cell is not unique and must be dynamically reconstructed 
throughout the simulation. In CFD-ACE-i-, this is accomplished through an upwind scheme with 
PLIC (Piecewise Linear Interface Construction method— i^), taking into account surface tension 
to determine accurately the interface curvature. As to boundary conditions, a zero static pressure 
is imposed at both inlet (liquid only) and outlet (gas only). At the initial time both fluids are at rest, 
with a short portion of the channel at the inlet filled with liquid. This specific configuration allows 
the liquid- vapour interface to assume the curvature corresponding to the hydrophilic partial wetting 
condition imposed at walls through the contact angle 0, thus generating the capillary pressure that 
drives the fluid motion. The CFD simulation set-up consists of a 2D straight channel with a height 
H = 800 nm and overall length L = 80 jim. At the bottom wall, a squared post, of side and height 
h = 400 nm, is located. The computational domain has been discretized with a squared structured 
non-uniform grid consisting of 185000 cells and 190000 nodes. 

Lattice Boltzmann method for multiphase flows 

The Lattice Boltzmann (LB) method is based on a minimal form of the Boltzmann kinetic equation 
describing the time evolution of discrete particle distribution function fi{x,t), denoting the proba- 



(O = 0(0/ + (1 -^)COg 



(5) 
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bility of finding a particle at lattice site x and time t moving along the lattice vector c,- ( Figure 2 ), 
where the index / labels the discrete directions of motion. 




Figure 2: The two-dimensional, nine-speed LB scheme. 



In mathematical terms the LB equation reads as follows:— 



Mx + cAt,t + At)-Mx,t) = (^fi{x,t) ~f^"'\p,pu)^ +F, At 



(6) 



where p is the fluid density, u is the fluid velocity and z = — 8 labels the the nine-speed, two- 
dimensional 2DQ9 model. The second term on the right hand side is a model collision operator 
Q.i, describing the relaxation towards a local equilibrium on a time scale T. Finally, Fi describes 
the effect of external/internal forces on the discrete particle distribution. The macroscopic density 
and momentum are obtained as weighted averages of discrete distributions:— 



(V) 



For the case of a two-phase fluid, the interparticle force reads as follows: 



-'^c^Y^WiXir{x,t)\i/{x + CiAt,t)ci z = 0-8. 



(8) 



Here V^(jc,?) = i//(p(j:,?)) = (1 — exp(— p(jc, ?))) is the pseudo-potential functional, describing 
the fluid-fluid interactions triggered by inhomogeneities of the density profile, and ^ is the strength 
of the coupling (see^>^ for details). Finally, Fi = Wi F ■ Ci/cl, where w,- = 0, 1 /9, 1 /36 are the 
standard weights of the nine-speed two-dimensional lattice 2DQ9^^ and Cs = 1/v^ is the sound 
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speed in lattice units. Besides introducing a non-ideal excess pressure p* = ^^c^i//^, the model 

also provides a surface tension 7~ '^c'l ^^J^ , where 51// is the drop of the pseudo-potential across 

the interface of width 5^. By tuning the value of the pseudo-potential at the wall V^(Pw), this 

approach allows to define a static contact angle, spanning the full range of values G [0" : 1 80"^] .— 

As to the boundary conditions, the standard bounce-back rule is imposed: any flux of particles 

that hits a boundary simply reverses its velocity so that the average velocity at the boundary is 

automatically zero. This rule can be shown to yield no- slip boundary conditions up to second 

order in the Knudsen number in the hydrodynamic limit of single phase flows. 

^ 2L ^ 

Effective channel 

(b) 

feriadic boundary conditions 

f -► Liquid; Pi Gas B Pg ^ [h 

t. ^ 

Figure 3: Geometrical set-up of the LB simulations. 



In this work we consider a ID channel composed of two parallel plates separated by a distance 
H = 40A, where A is the space discretization unit. This two-dimensional geometry, with length 



2L = 3000A, is divided in two regions, as shown in Figure 3 The left part has top and bottom 
periodic boundary conditions, so as to support a perfectly flat gas-liquid interface, mimicking a 
"infinite reservoir". The actual capillary channel resides on the right half, of length L. The top 
and bottom boundary conditions are those of a solid wall, with a given contact angle 9. Periodic 
boundary conditions are also imposed at the west and east sides. The squared obstacle, of side 
h = H /2 = 20 A, is placed on the lower wall. 



Molecular Dynamics 

During the last few decades Molecular Dynamics has proved to be one of the most broadly used 
simulation techniques, capable of reproducing many details of macroscopic fluid dynamics. It has 
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been extensively used to study dynamic wetting problems and to shed more light on the molecular 
processes close to the contact line region.— An appealing feature of MD (e.g., compared with 
Monte Carlo or molecular statics) is that it follows the actual dynamical evolution of the system. 

We perform MD simulation of a simple generic model on a coarse-grained level. In our MD 
simulation, see Fig. 4, the fluid particles interact with each other through a Lennard- Jones (LJ) 
potential. 



where e = 1.4 and a = 1.0. The capillary walls are represented by particles forming a triangular 
lattice with spacing 1 .0 in units of the liquid atom diameter o. The wall atoms may fluctuate around 
their equilibrium positions, subject to a finitely extensible non-linear elastic (FENE) potential 



withi?o =1-5 and = LOksT where ks denotes the Boltzmann constant and T is the temperature 
of the system; r is the distance between the particle and the virtual point which represents its 
equilibrium position in the wall structure. The FENE-potential acts like an elastic string between 
the wall-particles and their equilibrium positions in the lattice and keeps the structure of the wall 
densely-packed hexagonal. In addition, the wall particles interact with each other via a LJ potential 
with Cvvw = 1.0 and o^w = 0.8. This choice of interactions guarantees no penetration of liquid 
through the wall, and at the same time the mobility of the wall particles corresponds to the system 
temperature. 

Molecules are advanced in time via the velocity-Verlet algorithm with integration time step 
5t = O.OUq where to = (a^m/48e/j) = 1/ -\/48 is the basic time-unit and we have taken parti- 
cle mass m = 1 and kgT = 1. The temperature is maintained by a Dissipative-Particle-Dynamics 
(DPD) thermostat, with friction parameter ^ = 0.5, thermostat cutoff rc = 2.5a and step-function- 
like weight functions.— 1^ Fluid properties and flow boundary conditions arise then as a conse- 
quence of the collision dynamics and the local friction controlled by the DPD thermostat. An 



C/i,(r)=4e[(a/r)i2-(a/r)6] 



(9) 



(10) 
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advantage of the DPD method in comparison with other MD schemes is that the local momentum 
is conserved, so that the hydrodynamic behavior of the liquid at large scales is correctly repro- 
duced. 

It is worth emphasizing that all contact angles 6, used in the simulations, have been determined 
by measuring the mean static contact angle of a sessile meniscus in a slit, formed by two parallel 
atomistic walls for a few given liquid -wall interactions £. Specific contact angles are then chosen 
by interpolation between the relevant values of e. Therefore, the actual value of the contact angle 
when the fluid front is located at an obstacle edge, for instance, may incidentally deviate from 
the nominal value of . A recent MD study— has shown that the LW-law (Eq.(l)), holds almost 
quantitatively down to nanoscale tube diameters. In this study we use an obstacle shaped as a 
rectangular ridge which runs perpendicular to the flow direction in the slit and has the same atomic 
composition as that of the slit walls. The height and width of the ridge are chosen as lOcr so that 
the obstacle height takes half of the slit thickness. 



Numerical results 

We have performed simulations of spontaneous capillary filling in presence of a squared obstacle, 
for each of the three computational methods described above. The test case is the same: the filling 
of a channel of height H with a squared post of side h = H/2. In LB and MD, the post is placed at 



the central position on the bottom wall( Figure 1 1, at a distance / = 25H and I = 5H from the inlet. 



respectively, due to computational costs. Since the LW equation does not depend on the channel 
length, this is not a limiting factor for the present study (see below). Since the inertial transient 



with the CFD method is much longer than the estimated inertial time, as shown also in Figure 5 
in this case the channel length has been taken / ^ 100//, so as to guarantee that the front meets the 
post when the asymptotic regime is attained. 

We have performed a series of four simulations (for each method), varying the contact angle 
6 from 30° to 64°. The specific values of the physical parameters for each simulation method are 
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reported in Table 1 . 

Table 1: Physical quantities for the three different cases. Units have been chosen as follows: 
LB) Ax = H/NY, At = Vlb Ax^/v, Am = pi AjcVpls); MD) o = H/20, x = Vmd o'^/v, 5m = 
pi (7^/pMD- The different values of the gas density and surface tension reflect the computational 
constraints of the different methods. 



Parameter 


CFD value 


LB value 


MD value 


Channel height H (m) 


8-10-' 


8 -10-' 


8- 10-' 


Water density pi (kg/m-^) 


10^ 


10^ 


10^ 


Water kinematic viscosity v (m^/s) 


10^ 


10^ 


10^ 


Water dynamic viscosity rj (kg/(m s)) 


10-3 


10-3 


10-3 


Liquid/gas surface tension y (N/m) 


7.2-10-2 


9.6- 10 2 


3.4-10-4 


Vapour density pg (kg/m^) 


1.167 


28.12 


10.6 



Some comments are in order. By construction, the CFD approach is able to reproduce exactly 
all the physical properties of the flow. On the other hand, both LB and MD show some discrep- 
ancies in the value of the surface tension and vapour density. However, these discrepancies have 
been found to be irrelevant to the purpose of investigating the macroscopic features of capillary 
imbibition'^^^ 

As is well known, depending on the value of the contact angle 9, two scenarios may appear 
according to the Gibbs, or Concus Finn criterion:— "^^^ a) for small contact angles the front is able 
to climb the obstacle, walk on it, and eventually pass it; b) for large contact angles the front climbs 
the obstacle, but pins at its back edge, thus stopping the fluid motion. We wish to emphasize that 
the Concus-Finn, or Gibbs, criterion is based upon thermodynamic arguments and, consequently, 
it has been rigourously proven only at a macroscopic level. 

Recovering the Lucas- Washburn regime 

As is well known, the Lucas-Washbum asymptotic regime sets in once inertial effects die out and 
steady propagation results from the sole balance between capillary drive and dissipation. During 
the transient regime, the dynamic contact angle, which depends on the front speed itself, changes 
in time until its static value is attained.-^^*^ 



12 



Panillarv fillino wrfVi wall rnmioafinns 



100^ A 



10 



z/H 



CFD 30 
LB 30 
MD 30 
LW30 




100^ 



10 



CFD 40 
LB 40 
MD 40 
LW40 





Q J I I I I I I g ^ I I I I I I 

1 10 100 1000 10000 le+05 ' 1 10 100 1000 10000 le+05 

100 



10 



z/H 




0.1 



I 



10 100 



® CFD 50 

H LB 50 

A MD 50 

-■■ LW50 




0.1 



I I 



Jo le+05 1 



t/t 



cap 



10 100 
t/t 



CFD 60 
LB 60 
MD 60 
-- LW60 







)0 le+05 



cap 



Figure 4: Log-Log plot of the dimensionless front coordinate z{t) /H \s the dimensionless time 
t /tcap- Here z denotes the centerline position (y = H /2) of the front, y being the cross-flow coor- 
dinate. All simulations show superposition with the LW-prediction before reaching the obstacle, 
at all inspected contact angles. The arrows indicate the pinning points for the case of 50° and 60°, 



not visible on the scale of this figure. See also Figure 5 , Figure 6 and Figure 7 for a close-up of the 
dynamics around the obstacle region only. 



In Figure 4 the time evolution of the centerline (y = H/2) front position is reported as a func- 
tion of time, and compared with the dimensionless LW-law Eq. As one can see, all three 
methods exhibit satisfactory agreement with the theoretical LW prediction, although on a different 
range of time-scales. This is due to the different values of the parameters used in each method, and 
particularly, to the fact that the MD capillary speed Vcap = 7/^7 is 20 -i- 30) times higher than in the 
other methods. Indeed, it should be noted that the LW asymptotic solution sets in after a typical 
transient time of the order of T ~ //^/12v/, or in units of capillary times, T/tc-a„ = Based on 



the data in Table I, one can readily check that x/tcap ~ 1 for LB and CFD, and x/tcap ~ 10 ^ for 



MD. It is therefore possible to appreciate the anomalous transient in the CFD case, see Figure 4 
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which is of the order x/tcap ~ 100. However, since CFD and LB simulations last over 10^ cap- 
illary times, and MD simulations last about 10 capillary times, before reaching the obstacle, the 
superposition with the LW regime is free from transient phenomena. 



Front morphology while crossing the obstacle 



In Figure 5[ Figure 6 and Figure 7 , snapshots of the fluid front during the surmounting of the 
obstacle are reported for the three simulation methods (the figures report density contours). As one 
can see, the front dynamics and morphology are strongly affected by the presence of the obstacle: 
the liquid impinging on the obstacle must adjust to a 90° degrees discontinuity of the contact angle, 
which is clearly causing a significant deformation of the front, before the static value is recovered 
again on the flat-top surface of the obstacle. These changes of shape are well visible in figures 
Figure 5 and Figure 6[ for both CFD and LB. The case of MD is less clear-cut, due to the absence 
of a well-defined interface and to molecular fluctuations. Indeed, although the fluid meniscus 

one 



is clearly visible, its surface appears rough and strongly fluctuating in time. From Figure 7 



can see that some atoms evaporate from the liquid and overcome the slit. Moreover, long before 
the fluid meniscus has passed the obstacle, vapor condensation and partial filling of the wedge. 



formed by the rear wall of the ridge and the slit wall, take place (see also Figure 7i). Since fluid 
imbibition in a capillary is accompanied by the faster motion of a precursor film far ahead of the 
meniscus, -^^^"^ it is also conceivable that the wedge is filled by this atomically thin precursor. 
Since the meniscus position at time t is difficult to locate precisely, we rather measure the volume 
of the fluid in the capillary which is readily obtained by the total number of fluid atoms, residing 
at time t in the slit during the filling process. For an incompressible fluid, this volume is directly 
proportional to the distance z{t), travelled by the fluid front at time t. The curve giving z{t) is 



shown in Figure 7 d for several values of the contact angle, 30° < 9 < 64°. The shape changes 
also entail a substantial change of the front speed, as well documented in Figure 5\ i, Figure 6 3 and 



Figure 7 d , from which the front speed can be read off simply as the slope of the curves reporting 
the front position as a function of time. These figures show evidence of a significant acceleration 
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of the front (centerline position y = H /2) m the climbing stage, followed by a deceleration in the 
stage where the front is approaching the rear corner, where the chance/risk of pinning is highest. 
Here the fate of the front becomes critical. According to the CF analysis, if the contact angle is 
below 45°, there is enough drive from the upper wall to pull the front away from the rear edge of 
the obstacle. Otherwise, the front stops moving. 

The CFD and LB simulations confirm this picture, showing evidence of pinning only for the 
two angles above 45°. More precisely, they show a significant front acceleration in the climbing 
stage, followed by a coasting period, once the rear edge is reached. For = 30° and = 40°, the 
coasting period exhibits a finite lifetime, after which the front regains its motion. For 6 = 50° and 
6 = 60°, the coasting period does not seem to come to an end (within the simulation time), and the 
front is pinned. After the obstacle, the (unpinned) front is slowed down, as one can appreciate from 
the slightly reduced slope of the front dynamics, see Figure 5\ ) (only for 6 = 40) and Figure 6 



The MD simulations, on the other hand, tell a different story, in that even at = 50°, the 
front proves capable of overcoming the obstacle, although with a drastically reduced speed. Note 
that, while both CFD and LB show evidence of a strong front acceleration as they approach the 
obstacle, MD simply shows a monotonic deceleration. This is due to the fact that while in CFD 
and LB simulations we measure the distance travelled by the interface midpoint, as already pointed 
out, MD measures the total volume of the fluid in the capillary. 

Given the qualitatively different outcome of CFD(LB) versus MD simulations at > 45°, we 
next discuss the dependence of the filling dynamics on the different contact angles, case by case. 

• Contact angles 30° and 40° 

As already mentioned, in these cases, all three methods give the same qualitative outcome: 
in proximity of the ridge, the front deforms in response to the geometrical discontinuity, 
climbs up the obstacle and walks over its top. Once the rear-edge is reached, the bottom 
end of the front pins at the comer, and keeps moving on the top wall. Then, according to 
the CF criterion, the front overcomes the obstacle and completely fills up the entire channel. 
Manifestly, the standard z versus t relationship, described by the LW-law © is strongly 
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Figure 5: a) Snapshots of the front dynamics at different stages of a post surmounting CFD simu- 
lation. The various columns from left to right correspond to contact angles = 30°, 40°, 50°, 60°, 
as indicated, b) Time evolution of the interface midpoint while crossing the obstacle in CFD sim- 
ulations. The inset shows the instantaneous propagation speed of the front for the case G = 30°. 



violated in the vicinity of the obstacle. After overcoming the obstacle, the usual LW regime 
is recovered, although after a transitional period of time, which depends on the wettability 
6,— and with a reduced velocity. In particular, in all cases, the front is slowed down right 
after passing the obstacle, see Figure 5p , Figure 6 3 and Figure 7 3. However, at 6 = 30, 
the asymptotic behaviour is quite different for the three methods: (i) In CFD simulations. 



Figure 5 3, after a short transient time, the velocity basically regains the initial value before 



passing the obstacle, (ii) LB shows a similar behaviour, with only a slight reduction of 
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Figure 6: a) Snapshots of the front dynamics at different stages of a post surmounting LB simula- 
tion. From left to right, the columns correspond to contact angles = 30°, 40°, 50°, 60°. b) Time 
evolution of the interface midpoint while crossing the obstacle in LB simulations. The inset shows 
the instantaneous propagation speed of the front for the case = 30°. 

the front velocity after passing the obstacle, (iii) in MD, however, the front undergoes a 
substantial velocity decrease. 

This might be due to the details of the fluid- solid interaction during the obstacle surmounting 
(in which the contact angle is bound to undergo drastic changes), as well as to the different 
time-scales. A detailed analysis of this effect shall be deferred to a future study. 



Contact angle 50° 

In this case, while results from CFD and LB confirm the CF criterion, MD simulations show 
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Figure 7: a) Variation of the number of fluid particles in the capillary with time t for contact 
angles 9 = 32°, 40°, 50°, 64°. Symbols represent MD results, whereas lines denote the predicted 
behavior according to the LW-equation. The vertical grey-shaded column indicates an extension in 
the time axis, b) Time evolution of the interface, computed through the fluid volume measurement, 
while crossing the obstacle in MD simulations. The inset shows the instantaneous propagation 
speed of the front for the cases 9 = 32° . 



a different behaviour for the front motion: interestingly, the Concus-Finn (Gibbs) criterion 
for contact line pinning at the edge of the ridge is found to break down. Indeed, for 9 = 50°, 
the fluid front overcomes the obstacle in manifest violation of the Concus-Finn criterion. In 



order to vividly explain this feature, we show in Figure 8 the snapshots and in video 1 and 2 
the movie of the front dynamics, respectively for the cases 9 = 40°, 50°. These observations 
suggest that, on the nanoscale, the overcoming of topological obstacles is strongly affected 
by interface fluctuations, thus undermining the deterministic nature of the imbibition pro- 
cess. Indeed, both the CFD and LB method work in absence of fluctuations, and this would 
explain the difference with MD simulations. The problem of contact line pinning during 
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capillary imbibition acquires thus a stochastic character and is most probably governed by 
the size of the obstacles around the CF critical point. In fact, depending on the height of 
the ridge obstacle, a coalescence of the pinned meniscus with the molecules ahead of the 
obstacle, in the vicinity of the edge, may occur at later times. 

• Contact angle 60° 

Again, in this case all the three methods give the same result: the front deforms, climbs the 
obstacle, walks on its top, but pins at the back edge and definitely stops moving. 

e = 40° e = 50° 

time time 




Figure 8: Snapshots from MD simulations with d = 40° (on the left) and for 6 = 50° (on 
the right). It is clearly seen that for both values of the contact angle the front overcomes the 
obstacle, although at different times. This shows that the Concus-Finn criterion is violated 
in the case with = 50°. 



Further discussion 

In order to gain a better understanding of the previous results, and particularly of the viola- 
tion of the CF criterion for mildly super-critical angles in MD simulations, additional runs 
have been performed. More precisely, we have run MD simulations at 6 = 50° with a taller 
obstacle, h = 15(7, in order to inspect whether the front is still capable of surmounting the 
obstacle in violation of the CF criterion (in this case, the total slit height was correspond- 
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ingly increased, so as to leave the same clearance above the obstacle as in the previous 
simulations). Moreover, we have also run additional CFD and MD simulations at 6 = 50° 
with shorter obstacles, h = H /A, in order to inspect whether even CFD(LB) would show 
violations of the CF criterion upon reducing the obstacle size. The main outcome is as fol- 
lows: MD simulations with h = 15a do show front pinning, indicating that violations of the 



CF criterion disappear once the obstacle is made sufficiently tall (see Fig. Figure 9). This 
corroborates our previous conjecture of a (Arrehnius-like?) dependence of the CF criterion 
on the obstacle height, in the presence of molecular fluctuations. At the same time, CFD and 
LB simulations with h = H /4at 6 = 50° keep showing evidence of pinning, thereby lending 
further support to the idea that the CF criterion remains insensitive to obstacle size, so long 
as molecular fluctuations can be neglected. 



(a) 




20 40 60 80 
t/L 



100 120 



Figure 9: Snapshots (a) and centerline front coordinate (b), from MD simulations with 
6 = 50° and h = \5o. Unlike the same case with h = lOo, the front is now pinned in 
accordance with the CF criterion. 



Before concluding, a few words of caution should be spent on the fact that, being diffuse- 
interface methods, both CFD and LB propagate a finite-width interface. As a result, one may 
wonder whether and to what extent finite-width effects could interfere with the physics of 
post surmounting. Indeed, as shown in previous studies by these and other authors (— i^), 
finite interface width eventually leads to mild deviations from the LW law. However, they 
do not affect the qualitative outcome of a pinning/no-pinning test, unless the interface width 
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becomes comparable to the characteristic obstacle width. Both CFD and LB simulations are 
visibly far from this critical limit, which is why we are confident that the qualitative con- 
clusions of the present work are not significantly affected by finite- width effects. One may 
wonder whether such deviations might be paralleled to the effect of molecular fluctuations. 
To this regard, it is worth underlining that such a parallel has to be taken very cautiously, 
since molecular fluctuations stem from the microscopic physics of the problem, while finite- 
width effects are due to a lack of numerical resolution. Sometimes a mapping between the 
two can be established, but this can by no means be taken as a general rule. Finally, another 
potential source of discrepancy between CFD(LB) and MD is the fact that the latter allows 
slip motion while CFD and LB (in this set up) do not. Although any solid statement on the 
effect of slip flow on the dynamics of post surmounting must necessarily be deferred to a 
detailed quantitative analysis, we observe that due to the weakness of inertial effects, hy- 
drodynamic boundary conditions have little or no effect, on the dynamics/energetics of the 
obstacle surmounting. 

This is confirmed by the fact that slip flow effects should manifest through visible deviations 
from the LW regime, of which we have no evidence, at least in the parameter regime investi- 
gated in this work. This situation can drastically change in the presence of superhydrophobic 
effects, although we shall not be concerned with these problems in the present work. 

Conclusions 

Summarizing, we have studied the effect of geometrical obstacles in microchannels on the process 
of capillary filling, by means of three distinct simulation methods - Computational Fluid Dynam- 
ics (CFD), Lattice-Boltzmann Equations (LBE) and Molecular Dynamics (MD). The numerical 
results of these approaches have been compared and tested against the Concus-Finn (CF) criterion, 
which predicts pinning of the contact line at rectangular ridges perpendicular to flow for contact 
angles 6 > 45°. While for 6 = 30° (flow), = 40° and 6 = 60° (no flow) all methods are found 
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to produce data consistent with the CF criterion, at 9 = 50° the numerical experiments provide 
different outcomes. While pinning of the liquid front is observed in both LB and CFD simulations, 
the MD simulations show that the moving meniscus overcomes the obstacle and the filling goes 
on, for a sufficiently small obstacle. This result indicates that the macroscopic picture underlying 
the CF criterion and hydrodynamic approach needs to be amended near the critical angle. Further- 
more, while in CFD and LB simulations the front re-emerges from the obstace surmounting with a 
nearly unchanged velocity, in the MD case the post- surmounting velocity appears considerably re- 
duced. These results suggest that, away from the critical value 6 = 45°, the issue of front-pinning 
in a corrugated channel can be quantitatively described by a kinetic Boltzmann approach or by the 
macroscopic CFD method. 

While the CFD software used in this work is well-suited to handle complex geometries, it also 
shows some physical and computational limitations, namely anomalous-long transients and the 
need of a large computational grid to assure the required accuracy, which entails a correspondingly 
long computational time, much closer to the MD requirements than to the LB ones. In the vicinity 
of the critical angle, the motion of the front exhibits a strong sensitivity to molecular fluctuations 
which cannot be accounted for by standard (non-fluctuating) LB methods, let alone continuum 
methods. In particular, the MD simulations show that molecular fluctuations allow front propaga- 
tion slightly above the critical value predicted by the deterministic CF criterion, thereby introduc- 
ing a sensitivity to the obstacle height (the CF criterion is restored for sufficiently tall obstacles). 
On the basis of the present results, it would be indeed of interest to explore whether fluctuating 
hydrodynamic methods, either in the form of stochastic hydrodynamics or fluctuating LB, would 
prove capable of reproducing the results of MD simulations.^^ Whether the probability of "tun- 
nelling" to the deterministically forbidden region 6 > 45" shows an Arrehnius-like dependence on 
the obstacle height, stands as an interesting topic for future research. 
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